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Synopsis 


In the theory of elliptic boundary value problems, the phenomenon called “fundamental 
solution” plays an important role. Fundamental solution and Green’s function are very 
closely related and, infact, they are almost one and the same. The difference between 
the two is that Green’s function is associated with a boundary value problem (i.e. with 
the differential operator and the boundary conditions), whereas fundamental solution is 
associated with only the differential operator. Once the Green’s function is constructed, 
the solution is readily represented in terms of the integral of the Green’s function and 
other datas of the problem. 

We know that the method of boundary integral equations is the starting point of 
the today’s boundary element method. Here, too, once the fundamental solution is 
introduced in place of the weighting function, we get a readymade representation of 
solution in terms of the boundary integrals. To that extent, the fundamental solution is 
very vital in the boundary element method. 

Our work in this thesis is mainly to provide a simple and effective numerical method 
for singularly perturbed elliptic boundary value problems using the boundary element 
method. This involves the construction of fundamental solution in a very special way 
so as the boundary element method suits to singular perturbation problems. 

In Chapter two, we start with Helmholtz equation involving a large parameter A. 
The fundamental solution of this Helmholtz operator is represented in terms of Hankel 
function. Since in the boundary element formulation, this large parameter appears 



V 


only in the fundamental solution, we go for an asymptotic expansion of the Hankel 
function, where we have taken only the leading term of the expansion. Using this 
approximate fundamental solution in the boundary element method we demonstrate 
that the singularly perturbed Helmholtz equation can be efficiently solved. Further, we 
have obtained the bounds for the error due the numerical integrations involved in the 
method and show that this bound tends to zero as A -4 oo. 

Chapter three, extends this idea to a rather more general linear singularly perturbed 
boundary value problems of elliptic type in two dimensions. This chapter is diveded 
mainly into two parts, where we discuss two different linear operators and show that, it 
is equivalent to considering the general linear problems of singular perturbation type. 
Here again we construct the fundamental solution for these operators. In this process, we 
make use of the already availabe fundamental solution of the Helmholtz operator. That 
is, this fundamental solution is extended in a classical procedure, which involves solving 
the corresponding Fredholm integral equation. This is done by successive approximation 
technique, and also shown that the resultant series is convergent for the large parameter 
A. An approximation procedure using the dual reciprocity method is presented to 
construct and evaluate the fundamental solution in an algorithmic way. Finally, we 
have performed some numerical experiments to show the efficiency of the method. 

In chapter four, we propose an iterative technique via linearization for non-linear 
singularly perturbed elliptic boundary value problems, wherein the linearized problems 
are solved by the approach discussed in the earlier chapters. The convergence of the it- 
erative technique is proved for very large A, and the computational experiments exhibits 
the power and simplicity of the method even for non-linear problems. 

Chapter five discusses an application of the method to a physical model in the 
Semiconductor Physics and which happens to fall under the category of the singular 
perturbation problems. That is, by proper scaling this non-linear coupled system is 



viewed as a singularly perturbed one. Again, via linearization, we decouple the system 
and solve it iteratively. The linearised problems are solved by boundary element method 
which uses the fundamental solution constructed in the earlier chapters. Computational 
results are comparable with already existing ones. 
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Chapter 1 
Introduction 


1.1 General Introduction of the Problem 


In this dissertation, we propose some numerical techniques for solving singularly per- 
turbed boundary value problems of elliptic type. We start with linear problems and the 
extend our work to non-linear ones. The second order elliptic linear singular perturba- 
tion problems can in general be described as follows. 
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It is assumed that the given data /, ai, 6i, Ci, 02, 62 and C2 are sufficiently smooth 
functions in Q 

The condition for the uniqueness of the solution of the above boundary value prob- 
lem Cl -f- Ac 2 < 0 is assumed to be satisfied in Q 

Wasow[65] writes that 

whenever in a differential equation the coefficients of the terms of the 
highest order are small by comparison with the coefficients of the other 
terms, the solution of the problem can be expected to show irregular- 
ities owing to the non-uniform dependence on the coefficients of the 
differential equation 


Let us indicate, in detail, why this irregularities happen. By setting e = 1/A, the 


boundary value problem (1.1)-(1.3) 

can be written as 



Lu = e^L^u Liu 

= 

in 0; 

0 < e < 1 

(1.4) 
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= u on 
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In order to solve this boundary value problem, one will be tempted to ignore the 
term “e^L2'u”, as e is very small. Then it will reduce (1.4) to 

Liu = f(x,y) in Q. (1.7) 

This equation will be called the ‘reduced equation’. Since it is a first order equation, 
the solution in general will not satisfy the boundary conditions (1.5)-(1.6) on whole F. 
That is, there will be a portion F* of F, where the solution of (1.7) will not satisfy the 
boundary condition, which means that in the neighborhood of F* in Q, the solution 
of (1.7) is not a good approximation. This forces us to think of some alternative way 
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to take some corrective measures in the neighborhood of F*. This involves the con- 
struction of the so-called ‘boundary layer terms’, which are asymptotically equivalent 
to zero everywhere in except for a small neighborhood of F*. Naturally, we call 
this neighborhood of F*, the ‘boundary layer’ region. This construction of asymptotic 
approximation of the ’boundary layer terms’ of the solution is usually the most difficult 
and most interesting part of the the analysis of singular perturbation problems. 

Singular Perturbation Probation problem can still be explained in a simple way. 
Consider a boundary value problem depending on a small positive parameter e. 
Under some conditions, one can construct the solution of by the method of per- 
turbation - i.e. as a power series in e with the first term being the solution of the reduced 
equation (1.7). For small values of e the series cannot often be presumed to converge 
uniformly in the entire domain Q. When such an expansion converges as e ^ 0, then 
the problem P^ is called regular perturbation problem. On the other hand, when does 
not have a uniform limit in as s -4- 0, then this straightforward perturbation method 
fails, and as a consequence one may miscalculate or even lose essential results. In this 
case, we call Pg a singular perturbation problem. Generally, when the small parameter 
is multiplied with the highest order terms of the differential equation, as in (1.4), then 
it happens to be singular perturbation problem. One can define a singular perturbation 
problem as a one in which no single asymptotic expansion is uniformly valid throughout 
the the domain as s 0. 

1.2 Asymptotic Analysis of Singular Perturbation prob- 
lems 

Singular Perturbation is now a maturing mathematical subject with a fairly long his- 
tory and strong promise for continued important applications throughout Science and 
Engineering. Though the basic intuitive ideas involving local patching of solutions can 
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be found in early work by Laplace, Kirchoflf and others, Prandtrs[55] paper at the 1904 
Leipzig mathematical Congress began the study of the fluid dynamical boundary layers 
by analyzing viscous incompressible flow past an object as the Reynolds number. Re, 
becomes infinite. 

The Navier-Stokes equations, in relation to the Reynolds number, have become one 
of the most striking examples of singular perturbations, leading to the idea of boundary 
layers, due to Prandtl. Right afterwards, many fluid dynamicists have got interested 
and started working on it. However, Friedrichs and Wasow seem to have been the 
first mathematicians to initiate the mathematical investigation of asymptotic solution of 
singularly perturbed boundary value problems.Their work was motivated by an analysis 
of the edge effect for buckled plates. They are the first one to use the term “Singular 
Perturbation” in the title of[24]. 

Soon afterwards, Levinson[45] began the study of wide spectrum of important top- 
ics in singular perturbations and made intuitive contribution. Mainly, his work on 
asymptotic properties was based on a maximum principal that can be derived for the 
differential equation (1.1). Later, Visik and Lysternik[64] simplified Levinson’s work, 
based on norm estimates. 

Beginning around 1950, fluid dynamist solved many interesting physical problems, 
such as the Linoleum-rolling problem (Carrier[10]) and low-Reynolds-number flow past 
bodies. At Caltech’s Guggenheim Aeronautical Laboratory, Lagerstrom, Cole, Latta, 
Van Dyke[62], Kaplaun and others became equally involved in asymptotic expansion 
procedures for more general singular perturbation problems, 

An oversimplified matching procedure was presented in the book of Van Dyke[62]. 
The straightforward recipe he provided made it easy for tremendous variety of scientists 
to learn the rudiments of matching and apply it to important problems in their own 
disciplines. The basic idea, much as in Friedrichs’ earlier work, involved an asymptotic 
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match of the inner and outer expansions at the edge of the boundary-layer, where they 
should both be appropriate. 

Cole[ll] stressed limit process expansions and two timing in a context for broader 
than fluid mechanics. Indeed, the results obtained through matching generally coincided 
with those known through the interactive folklore of the various fields. Eckhaus and 
DeJager[20] have presented a more complete theory of singular perturbation problems 
for second order linear differential equations of elliptic type. 

By 1970, courses in perturbation methods became common in Science, Engineer- 
ing and Applied Mathematics departments, and inevitably a string of textbooks and 
high-level monographs began to appear. They include Nayfey[51], Van Dyke[62_. Bell- 
man[4], Eckhaus[18, 19], Eckhaus k DeJager[20], 0’Malley[52], Hemker[31], Kevorkian 
k Cole[41], Kaplaun[39], Meyer k Parter[46], Erdelyi[22], Ames[l], Hughes[21] 

1.3 Numerical Analysis of Singular Perturbation Prob- 
lems 

The advent of sophisticated computing machines has opened up new areas in almost 
every field and singular perturbation problem is no exception. 

When an asymptotic analysis is valid and a few terms in a asymptotic expansion 
describe the solution sufficiently accurately, one usually can rely on standard techniques 
to obtain the solutions. However, when an asymptotic analysis is difficult to handle or 
performs badly - and this can happen very often considering the complexities occur 
in the practical problems - one then asks for a numerical method to solve the singular 
perturbation problem. In a sense, numerical methods are intended for broad classes of 
problem and are intended to minimize demands upon the problem solver; where as, the 
asymptotic methods treat comparatively restricted class of problems and require the 
problem solver to have some understanding of the behavior of the solutions expected. 
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With the availability of computers, the interest in the area of numerical analysis of 
singular perturbation problems is increasing among the applied mathematicians. Nu- 
merical treatment of singular perturbation problem has always been far from trivial, 
because of the boundary-layer behavior of the solutions, i.e. if the solution of a perturb- 
ation problem exhibits rapid variations within some small region whose extent tends to 
zero as the small parameter e — > 0, then the solution is said to posses boundary layer 
behavior. 

It is indeed a big challenge to provide a numerical method to a singular perturbation 
problem and the numerical analysts have accepted it. This is very evident from the fact 
that there already exists a vast literature in the field of numerical analysis of singular 
perturbation problems (for instance, see the survey article by Kadalbajoo k, Reddy [38]). 
However, all these people have mostly confined themselves to only ordinary differential 
equations and not many research work is done on numerical analysis of singularly 
perturbed partial differential equations. Of course, we do find research articles here 
and there in this direction but, overall, we feel that it has received very little attention. 

In the conference on Numerical Analysis of singular perturbation problems held 
at the university of Nijmegen, The Netherlands, there are about six papers that talk 
about partial differential equations. The first one is Brandt[7]. The approach he has 
proposed is to first discretize the given singularly perturbed boundary value problem 
using the finite-element or finite-difference schemes and then use the multi grid adaptive 
technique. This process automatically takes care of boundary layers. It is shown that 
the global error decreases exponentially as a function of overall computational work, in 
a uniform rate independent of the magnitude (s) of the singular perturbation terms. 

Griffiths[29] considers 

^ = (x.t)e(0,l)x(0.1) and A/e » 1 


.(1-8) 
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Using the semi-discrete generalized Galerkin method the problem is converted into 
a system of ordinary differential equations. Analysis of this system has been carried 
out from the point of view of phase errors, decay rates and the presence of otherwise 
unphysical oscillations in the solutions. 

Veldhuizen[63] first develops a method to solve a first order ordinary differential 
equation using the finite element space consist of piecewise quadratics and then extend 
it to a particular type of singularly perturbed second order elliptic boundary value 
problem. 

Axelsson and Gustaffson[2] have presented an interactive method that is essentially 
based on approximate factorization of the matrix in combination with a minimal residual 
conjugate procedure to solve the resultant linear systems that is obtained from the 
discretization of a singularly perturbed second order elliptic boundary value problems. 

Baranger[3] has come with a very interesting method that numerically calculates 
the thickness of the boundary layer. His work was motivated from the result that the 
solution of the 2-D boundary value problem 

— eAue -f Ue = u in 
Ueir = 0 

has an asymptotic expansion of the form 

(1.9) 

where T is taken to be a unit square and d(a:,r) is the distance from x € to T. 
This leads to 

d{x.V) = ^log{cs)-^ = 0{e^lHog{E)-^) (1.10) 

But since it imposes strong assumptions on V (going from to G®!) they propose 
a method to define the thickness of boundary layer based on estimates outside the 


Ur 


u + exp 
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boundary layer which is valid for a Lipschitzian boundary for second order operators 
with L°° coefficients. 

Inayat Noor and Aslam Noor[34] prove that by an iterative scheme obtained as 
a projection of a real Hilbert space into a convex subset the solution of variational 
inequalities can be found. The asymptotic solution obtained by this iterative method 
does converge strongly to the exact solution. And it is shown that this method can be 
extended to a class of singularly perturbed boundary value problems by satisfying some 
extra conditions. 

Jordan[37] in one chapter of his Ph.D. thesis, propose a numerical method for an 
elliptic boundary value problem of the form 


-sA + — + c(x, y) 


u 


u 


-f in 


where = {(a;, y) € : \x — y\ < 1}. 

The method is based on the asymptotic behavior of the solution; i.e., the solution 
admits an asymptotic expansion of the form 


u, = U + V + Z 


where U is the solution of the reduced problem. V is the correction to [/ in order to 
meet the remaining conditions and Z is the remainder term that satisfies the estimate 
jlzll = 0 (e) uniformly in D,. 

Schatz and Wahlbin[59] investigate the behavior of the Galerkin finite element method 
without any modifications, on a singularly perturbed second order elliptic boundary 
value problem. Some non-linear elliptic problem also has been discussed and finally 
some numerical experiments have been carried out. 
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Ramanujan and Sunanda Kumari[56] propose the method of lines(MOL) for a non- 
linear elliptic boundary value problem. First, it is shown that the ordinary method 
of lines doesn’t give satisfactory results for singularly perturbed non-linear boundary 
value problem. Then the modified method is proposed and prove its convergence. 

Goldstein has published a string of papers[28, 25-27] on Singularly perturbed (ini- 
tial and) boundary value problems. His work essentially presents a estimation of the 
condition number k{L\~^L\) in terms of the large parameter K and Ki where L\ and 
Li are finite element discretization of two singularly perturbed elliptic operators or 
order one and two respectively and h is the discretization parameter. Shown how to 
construct the preconditioning operators such that the condition number of the systems 
is bounded as the number of unknowns increases and, at worst, grows slowly with K 
and Ki increase. Numerical experiments have been performed on elliptic and parabolic 
convection-diffusion models and two dimensional Helmholtz problems. 

Boglaev[6] has presented parallel iterative algorithms via domain decomposition for 
singularly perturbed nonlinear problems. The domain decomposition is done according 
to the singular perturbation character of the problems. 

1.4 On Boundary Element Method 

Integral equation techniques in boundary value problems have been around for a very 
long time. As far back as 1903, Fredholm used discretized integral equations in po- 
tential problems which formed the basis for the “indirect” boundary element approach. 
This approach is referred to as indirect because it uses fictitious density functions or 
sources that have no physical meaning but can be used to calculate physical quantities 
such as displacements and stress. An integral equation relating boundary values of 
displacements and tractions had been established by Somigliana in 1886. Somigligna’s 
identity forms the backbone of the “direct” boundary element formulation. 
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Several books and papers on integral equations in potential and elasticity theory 
have appeared regularly by distinguished mathematicians such as Kellog[40] in 1929, 
Muskhelishvili[50] in 1953, Mikhlin[47] in 1957 and Kupradze[42] 1965. The integral 
formulations, however, were solved by analytical procedures that limited their applic- 
ations to simple problems. Until the sixties, there was no major work that extended 
the range of problems that can be solved by integral equations nor considered devising 
numerical algorithms for solving integral equations. 

In the early sixties high speed digital computers and numerical techniques started 
to find their way into engineering applications. In particular, the finite element method 
attached a great deal of interest and demonstrated how a wide range of complex engin- 
eering problems can be solved using the computer with an impressive accuracy. 

The main breakthrough in boundary integral solutions came in 1963 when two clas- 
sic papers were published by Jaswon[35] and Symm[60]. Their approach consisted of 
discretizing the integral equations of two-dimensional potential problems governed by 
Laplace equation into straight line elements over which the potential functions are as- 
sumed constant over each element. The elements were described in terms of nodal 
points and the integrations performed using Simpson’s rule except for some singular 
integrals which were integrated analytically. Jaswon and Symm’s approach can be clas- 
sified as “semi-direct” because the functions used to formulate the problem are not 
fictitious and can be differentiated or integrated to calculate quantities. 

Other similar integral equation approaches were adopted by Jaswon and Ponter[36] 
for torsion problems involving shafts with different regular cross-sections. Hess and 
Smith[32] for potential flow problems around arbitrary shapes and Harrington et a/[30] 
for two-dimensional electrical engineering problems. 

The first paper to use the direct approach of using displacements and tractions 
in an integral equation applicable over the boundary was published by Rizzo in 1967. 
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Rizzo’s work was very original in that it was the first to exploit the strong analogy 
between potential theory and classical elasticity theory and devise a numerical approach 
of solving th problem. He used straight line elements to discretize the boundary where 
the functions are assumed constant over each element. Simpson’s rule was used for all 
but the singular integrals. 

The extension of the direct integral equation approach to three dimensional problems 
by Cruse[13] in 1969 followed a very similar formulation to Rizzo’s works except that 
the surface was discretized into flat triangular elements with the displacements and 
tractions assumed constant over each element. Cruse’s work was continued in two 
further publications to include practical three-dimensional applications and comparison 
with finite element method (Cruse[15] in 1973) and using more sophisticated elements 
by allowing the variables to change linearly over each element (Cruse[16] in 1974). 

During that early period of development (from 1967 to 1972), the integral equation 
formulations were extended to include non-homogeneous problems containing inclu- 
sions (Rizzo and Shippy[58] in 1968), anisotrophic materials (Cruse and Van Buren[17] 
in 1971 and Cruse[l4] in 1972). These early publications were crucial because they 
provided the firm foundation for further boundary element developments and demon- 
strated that the boundary element approach was a very powerful and accurate numerical 
technique that should be taken seriously. 

Many of the algorithms and numerical modelling methods used in finite element for- 
mulations also found their way into boundary element formulations, such as the concept 
of higher-order elements, by using quadratic shape functions. These elements were first 
used by Lachat[43] in 1975 and improved further by Lachat and Watson[44] in 1976. 
Other authors, such as Tan and Fenner[61], used isoparametric quadratic elements, 
where both geometry and variables are allowed to change quadratically over each ele- 
ment, and demonstrated the high resolution of stress obtained in three-dimensional 
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problems. 

Since the early seventies the boundary integral equation approach has continued to 
develop at a fast pace and has been extended to include a very wide range of continuum 
mechanics including advanced non-linear applications. Over the years, the integral 
formulations have been referred to as either the boundary^ integral equation method or 
the boundary element method. 

Comparison of the features of the boundary element method with finite element 
method can be summarized as follows. 

• Less date preparation time. This is a direct result of the reduction of dimension- 
ality by one. Thus the analyst’s time required for data preparation for a given 
problem should be greatly reduced. Furthermore, subsequent changes in meshes 
are made easier. This advantage is particularly important in problems where re- 
meshing is required. 


• Less computer time and storage. For the same level of accuracy, the boundary 
element method uses a lesser number of nodes and elements. Since the level of ap- 
proximation in the boundary element solution is confined to the only boundaries, 
boundary element meshes should not be compared to the finite element meshes 
with the internal points removed. 

• Less unwanted information. In most engineering problems, the ‘worst’ situations 
(such as fracture, stress concentration and thermal shock) usually occur on the 
surface. In many design codes and engineering practices, the analyst is usually 
only concerned with what happens in the worst situation. Thus modelling an 
entire three-dimensional complex body with finite element and calculating stress 
at every nodal point is very inefficient because only a few of these stress values 
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will be incorporated in the design analysis. Therefore, using boundary elements 
is a much more efficient use of computing resources. Furthermore, since internal 
points in boundary element solutions are optional, the user can focus on a partic- 
ular interior region rather than the whole interior. 

More features of the Boundary Element Method can be found in Brebbia[8, 9]. 

Of course, there appears to be some disadvantages in boundary element method, one 
of which is that the solution matrix resulting from the boundary element formulations 
in unsymmetric and fully populated with non-zero coefficients. However, this is not a 
serious disadvantage because to obtain the same level of accuracy as the finite element 
solutions, the boundary element method needs only a relatively modest number of nodes 
and elements. 

1.5 Summary of the thesis 

Our work in this thesis is mainly to provide a simple and effective numerical method 
for singularly perturbed elliptic boundary value problems using the boundary element 
method. This involves the construction of fundamental solution in a very special way 
so as the boundary element method suits to singular perturbation problems. 

In Chapter two, we start with a singularly perturbed Helmholtz equation. That is, 
the equation involving a large parameter A. The fundamental solution of this Helmholtz 
operator is represented in terms of Hankel function. Since A is very large, we go for an 
asymptotic expansion of the Hankel function, where we have taken only the principle 
term in the expansion. Using this approximate fundamental solution in the bound- 
ary element method we demonstrate that the approach is very effective for singularly 
perturbed boundary value problems. 
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Chapter three, extends this idea to a rather more general linear singularly perturbed 
boundary value problems of elliptic type in two dimensions. This chapter is divided 
mainly into two parts, where we discuss two different forms of the equation (1.1). In 
the construction of the fundamental solution for these linear operators, we make use 
of the already available fundamental solution of the Helmholtz operator. That is, this 
fundamental solution is extended in a classical procedure, which involves solving the 
corresponding Fredholm integral equation. This we do by successive approximation 
technique, and show that the resultant series is convergent despite the presence of the 
large parameter A. An approximate procedure is presented to construct and evaluate the 
fundamental solution in an algorithmic way and subsequently performed some numerical 
experiments to show the efficiency of the method. 

In chapter four, we propose an iterative technique for non-linear singularly perturbed 
elliptic boundary value problems, wherein the linearized problems are solved by the 
approach given in chapter three. The convergence of the iterative technique is proved 
for very large As, and the computational experiments exhibits the power and simplicity 
of the method even for non-linear problems. 

Chapter five discusses an application of the method to the problems in the Semicon- 
ductor Physics. Here we deal with a singularly perturbed non-linear systems and the 
iterative technique proposed in chapter four is extended to this systems. Computational 
results are comparable with already existing ones. 



Chapter 2 

The Method for Singularly Perturbed 
Helmholtz equation 


2.1 Introduction 

To begin with, we take-up the Helmholtz equation involving a large parameter and 
propose to apply our technique to this boundary value problem. As we said earlier, the 
role of fundamental solution is very important in the boundary element method. Keeping 
the large parameter in mind, we take the asymptotic expansion of the fundamental 
solution of the Helmholtz operator and use it in the boundary element method. Error 
bounds of the numerical integrations involved in the method have been obtained along 
with the computational results. 

2.2 The Problem 
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where the boundary T = Fi + r 2 and u, q are known functions. 


Or, equivalently, by setting A = l/e, the problem can be restated as 


Au - }?u = 0 in n l<CA<co (2.4) 

u = u on Fi (2.5) 

« (= 1 ^) = ’ “ 

which is the well-known Helmholtz equation. 

This is a particular case of the equation (1.1). i.e., we have taken = A and 
Lj = -/, where I is the identity operator. 


2.3 Fundamental Solution 

As we said earlier, the concept of fundamental solutions is closely related with the 
Green’s functions. For a given differential equation 

Lu = f (2.7) 

suppose that L is invertible, i.e. there exist an operator L~^ such that L = LL ^ = 
/, where / is the identity operator. As L is a general differential operator, it is natural 
to expect to be an integral operator involving a kernel P{z,^), where z and ^ are 
points in Q.. Thus we write 

u{z) = = I -P(z, (2-8) 

Without loss of generality, we have introduced a negative sign in the above integral just 
for some convenience. Suppose now we operate oh both sides of (2.8) with L again. 
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Then 

L{u{z)] = f(z) = Lj -P(z, (2.9) 

As T is a differential operator depending on z, we can write L inside the integral. 

j LP(z,e)f(m=-n^) ( 2 . 10 ) 

which requires that 

iP(z,5) = -i(|z-® (2.11) 

where 5 is the Dirac delta function and I 2 — is the distance between 2 and And 
we call this the fundamental solution of the operator L. That is, if we replace 

/ by in (2.7) then the solution is called the fundamental solution. 

In case of Helmholtz operator, the fundamental solution is 

P(z,e) = (i/4)ifi^^(iA|z-e|) (2.12) 

where is the Hankel function of the first kind of order zero. 

2.4 Formulation of the Method 

In this section we present the boundary element method in its simplest form, i.e., the 
(boundary) elements are linear and the shape functions are constant over each element[8, 
9]. 

We integrate equation (2.4) over the domain with a weighting function w and get 

f [Aw — X‘^u]wdQ, = 0 
Jo. 

[ Au.wdQ — [ }?u.wdQ. = 0 
Jo Jo 



18 


By using the Green’s second identity 



dUm 

dw 

1 [w.Au — u,Aw]dCl = 

dn 



A»..da + 1 - 1 n.|^dr - ^ = o 

- lu.^^ = - (2.13) 

The weighting function w is chosen as the fundamental solution, i.e. the solution of 


[Aw - X^w] = -5z 

where 6- = - '?!) is the Dirac delta function and z is the fixed point where the 

solution is required. Note that we treat ^ = (^ 1 .^ 2 ) as a variable in fi. If r is the 
distance between z and then we have 

,^(r) = ^Hi'\iXr) 


where i = v^- Since A is very large we use the asymptotic expansion of the Hankel 
function [49]. 


HS'\iXr) cs e' 


'i{iXi — 7r/4)] 


LTfiAr 


( 1 / 2 ) 


where we have taken only the leading term of the expansion. This leads to 


w{r) = '-BQ\iXr) 


2^g(-Ar-i:r/4) 


iirXri 


( 1 / 2 ) 
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= _%(-»x/4) ^(-Xr) 

4 2 (^/ 2 ) v/ttXt 


g(-Ar) 

2.\/27rAr 


and 


^ _ -1 (H-2Ar)e-^'’ 

4\/27rA 

Now, equation (2.13) becomes, 



-u(2) 


/. 

L 


dw 

u—dT 

on 

dw 

li— dr 
an 



w 


w 


dUm 

dn 

dU-m 

dn 


dr 

dr 


(By the property of Dirac delta ) 


u{z) + / 


JT. f 

u-^dr — + 

dn Jr 


dUm 
w— 
r an 


dr 


(2.14) 


This expression is valid for any arbitrary interior point z of 0,. Since we do not know 
u on ra and [du/dn] on ri, we seek an expression for u at a boundary point. Let 
Cj be a boundary point of r. Let be a semi-circle centered at Zi and of radius p. L' 
b(‘ the remaining part of L such that 


iimr' = r 


Considering Lp -b L' as the boundary, the point Zi becomes an interior point. So, we 
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apply equation (2.14) for this point and take the limit as /? — >• 0. 


lim 


^i^z) + / 

= lim 

f dUm ^-p 

/ w — — aF 

Jr'+Vp dn 

p-*0 

yr'+Fp an 


H-lim/ u^dV + lim [ u^dr = lim [ w^dF + limf w^dT 
p~>o Jr/ an p~>o Jr^ an p-^o Jr' dn p^o Jr^ dn 


u{^i) + ^ 


dw 


dw 


u^dT + lim f u~dr 
r on p->o yfp on 


= [ w^^dr + lim f w^^dT (2.15) 
Jr dn P->o Jr„ dn ^ ^ 


p-^o Jr^ dn 


(2.16) 


Consider the second integral term on the left hand side of (2.15). 

dvj 


lim / u^dT = lim 
p->-o JVp dn p-^o 


dT 


.Applying mean value theorem for integrals, 


lim/' 


, u I -;:r- 1 dr = limu* { ^ | vrp 
P-^oJr, \drj^^^ P^o \dr j 


r-p 


,-(l + 2Ap)e(-^^) 

lim u 7 == — TTp = 0 

P-+0 4\/27rA 


where u* is the value of u at some point of Vp. Similarly, we have 


lim [ (w) 

P-^oJrr 


du 

2V^pV2 


Trp = 0 
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Hence equation (2.15) becomes, 




9Um 

It;— —dr 

on 


( 2 . 17 ) 


for any boundary point Zi. 


Now, naming the elements of the boundary as Fi, r 2 , ....rn,we have F = U”=i Tj. 

Let the value of u on Fj is Uj and that of [du/dn) on F^ is q^. Let be the mid-point 
of the element Ff Then u{z^) = u^. 


So, (2.16) can be written as 


Ui -I- 


r dw r 

/ u^dT = j 

./Uj-Fj on Ju 


dUm 
w— 
UjTj an 


dT 


Ui + 



dw 

an 



w 


dUfa 

dn 


dF 


Ui -I- 




( 2 . 18 ) 


Call Gij = /r^ wdT, Hij = Jr. {dw/dn)dT, and set 


Hij = Hij for i^j 

= 1-1- Hij for i = j 
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Then (2.17) can be written as 

= {2.19) 

J=l J=1 

Here we know some of Wj’s and some of g^’s. i.e. we have exactly n known quantities and 
n unknowns. Writing the expression (2.18) for i=l(l)n, (i.e.) apply equation (2.16) to 
every mid-point of the elements Fi, Fa, ....Fn we have the equations in the matrix-vector 
form 

HU = GQ 


where 


H = [H,,Un G = 

Q = [qi,q 2 ,--Qnf 


Rearranging this system with known quantities on one side and the unknowns on the 
other we get the linear system of equations 


Ax = F 


VV'e solve this system by Gaussian Elimination method. Now, we know all the uj s and 
q/s (i.e.) we know both u and (du/dn) on F. 


Now, for an interior point z we have equation (2.14) 

dw r dUn 
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u{z) = - / u^dT + / w^dT 
Jv on Jt on 


dur, 


u{z) = -/ 
Ju 


dw f 

u-^dT + / w — — (ir 


dn 


lujVj dn 


3=1 

We can thus find the solution at any interior point of Q. 


2.5 Error Bounds 


In equation (2.17) we do numerical integration to compute the quantities Gij&cHij. We 
try to obtain some bounds for the truncation error as a function of h, the length of the 
element, and K where we use K-point Gaussian-quadrature formula. First we define 
the quantities as the distances from z„ the mid-point of F,, which is the point 

under consideration, to the starting point and the end point of Fj (see Fig-2.1). 


Take the transformation, 

r = l[(rj|^ - -f (r^f + r|]^)] = ^ (at + p) 
where, a = r\f - and P = -f- 



Now, we have, 
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Figure 2.1; and r^f are defined 


= [ w{r(t))dr{t) 

= »(') I -i‘ 


Of 

/=i 


ir 


where the error term 


^ 2 ^'^-H 2 K ) i (-“‘Si 


We have 


u;(r) 


,(-Ar) 


2\/27rAr(^/^) 


(2.20) 


( 2 . 21 ) 



( 2 . 22 ) 


w{t) 


V2e-^^ {at + 
2 ^/^ 


We know that, 


fP 2 K) 


(t) 


V2 ^ /2K\ 

2V2ttX ^ \ / / 





{at + P)^ 


\/2 ^f 2 K\f-\' 


2 ^/^ 

X 


/=0 


/ _x l“£db£l 
a^e ^ 2 


-l-2K + l + \^ (ai + /j)(-i/2-2ic+o 


< a2A-g[-Ai£^] ^ 2ir\ f-X\ 

~ y^{at + /?)(V2) ^ w y V 2 j 


( ( 1 / 2 ) ^ 
\ + P)) 


y^{at + pyi-^) V 2 2{at + p)) 


Q, 2 J'(rg[-AW^] ^ 

y^{at + pYm)+’^K22K 


Max 

-1 < i < 1 



< 


( \ 
V22^\/4^y 


r Max p-A(2^)] 

(1 + A(a* + /?)]“■ 


■ Min 

-i<t<i 

{at + 
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Let M — (oit + /3) and m — Since, ai + /5 = 2r > 0 always, 

we have m > 0 and M > 0. Therefore, we get 


Max 

-1 < t < 1 



< 


_wz_^ 

22^^/4^; 


e-xm (1 + 

^2iC+(l/2) 


Suppose h is the (boundary element) discretization parameter, then a = < h, 

which implies that 


Max 

■1< t < 1 




f \ e-^’” (1 + \M^^) 

^ V22^\/4^j m2^+(V2) 


Substituting this in (2.19), we get 


\Eij\ < 


22A'-2(2ii:)! 


g-Am (1 ^ XM^^) 
^2A+{l/2) 


l-^bl < 2‘^^-\2K)\ V22t^V4^J m2t^+(i/2) 


We not(i that, 

1. I^ijl — > 0 as A -4 oo for any K and h. 

2. \Eij\ — )■ 0 as iT -4 cx) for any A and h. 

3. \Eyj\ — >• 0 as /i -4 0 for any A and K. 

Following the similar lines, we can show that the truncation error in evaluating Hij 


tends to zero too. 
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2.6 Test Example 


—e'^Au + u 


0; 0 < a:, y < 1 and 0 < e <C 1 


u = e + e the boundary 


We solve the above problem[59] whose exact solution is the boundary condition itself. 
W e hav(^ tested the method with 4, 8 and 16 elements of the boundary and to evaluate 
tli(' matrix elements Hi^ & Gy we have used 4-point Gaussian quadrature formula, 
(.laus.sian Elimination method (with partial pivoting) is used to solve the resulting linear 
system. 

In table-2.1, we have tabulated the sup norm of the error function for various f’s. The 
much desired aspect here is that the error function decreases with s. This is because, 
in ii.sirsg the <isymptotic form for the fundamental solution, though we take only the 
first, term of the expansion, the truncation error is of order e which tends to zero as 
j 0. Also we note that the solution itself converges to zero as e 0 which coincides 
with the solution of the reduced problem. 

2.7 Conclusion 

In this chiipter, we have shown that the boundary element method, with the suitable 
imxlihcations, ctm be used to solve singularly perturbed boundary value problems where 
tb'. etpiation involved, in this case, is Helmholtz equation. We have also obtained the 
(‘rror bounds of the numerical integrations involved in the boundary element method 
and shown that they tend to zero as the discretization parameter tends to zero for A 
large. It is worthy to note that we have used only the simplest boundary element 
method in the numerical example and can be very well improved by using some higher 
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Table 2.1: Sup norm of the error function 




No. of Boundary Elements 


A 

4 

8 

16 

32 

10-^ 

.4632x10-2 

.3994x10-2 

.1256x10-2 

.1113x10-2 

10-‘ 

.8241x10“^ 

.8239x10-^ 

.8230x10-^ 

.8226x10-^ 

1Q-.3 

.7411x10-20 

.7411x10-2° 

.7411x10-20 

.7411x10-2° 

io-‘ 

.1931x10-^°® 

.1931x10-1°® 

.1931x10-1°® 

.1931x10-1°® 


order eh'uu'nts. The method is very simple and does not need any knowledge about the 


heliavior of the solution. 



Chapter 3 


Extension of the Method to General 
Linear Problems 

3.1 Introduction 

Our aim in this chapter is to present a method for solving a singularly perturbed gen- 
I'ral linear (dliptic problem (1.1) by extending the technique used in the second chapter. 
First, w(^ shall discuss two different forms of equation (1.1) and then show that the 
apf>roach can he used to any general linear elliptic problems of the form (1.1). Again, 
th(^ approach will be to construct the fundamental solution in a suitable way and show 
(he efh'ctiveness of the method for singularly perturbed problems of general linear type. 
Throughout this chapter, a,b,c and / are assumed to be smooth functions in 0. We 
show that the ((onstruction of the fundamental solution for the operator L in (1.1), can 
be done by considering the following problems. 

Probleiir 1 


t!l 

d 0 

A + a-r- + 6— + c 
dz oy 

u — = 

= / 

in 


1 4; A < oo (3.1) 



u = 

= u 

in 

Ti 

(3.2) 


A 


= g 

in 

La 

(3.3) 
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We start with the fundamental solution of the Helmholtz operator, [A ~ A^], and then 
extend it to construct the fundamental solution of The process involves solving the 
( onesponding the Fredholm integral equation which is done by successive approxima- 
tion method. Ihe resultant series is shown to be convergent for large parameter A. 


Problem 2 


ru ~s\u + A“ 


a— +b— + c 
ox ay 


u = A^/ in Q,; 1 C A < oo (3.4) 


w 


rith boundary conditions (3.2) and (3.3). 


the construction procedure is same as in the case of problem 1. The only 
(iifFer(>nc(^ is that the kernel involved in the Fredholm integral equation is different. In 
thi.s f.ast' also, we have shown that the resultant series in solving the Fredholm integral 
ecjuation convergCvS uniformly in Q. for large A. 

.Sectitui 3.2 show.s that the fundamental solution of L' can be done from [A — A^] 
and section .3.3 gives the fundamental solution for L". Combining these two, we can see 
that rite fimdamental solution can be obtained for the operator L from the fundamental 
solution of L" t)y following the lines of section (3.2). 

3.2 Construction of the Fundamental Solution for L' 

In this sectioti we coustruct[12, 57] the fundamental solution for the operator 


L'~ 


A + Cl 


du du 


+ c 


-A2 
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The fundamental solution of the operator [A - A^] is 

Pojz - e; A) ^ ^ [5_P,l-,lh- lil] ( 3 . 5 ) 

2v^ [ J 

From (3.5) we shall extend to the fundamental solution for the operator L'. In order to 
do this, let us seek a fundamental solution with a singularity at a point ^ = (^i, ^ 2 ) € Q 
in the form 

P(;,?,A) = P„(z-e.f.A)+ f P„(z->).i),A)A(i 7 .f,A)<i,n (3.6) 

J m 

where s, A) is determined such that when z ^ P{z,^,X) satisfies L'P = 0 . 

'Fhen from (3.6), 

VP„(z-^,i,\)+L' j Pa(z-n, V, A)A(' 7 . f , A)c(,n = 0 (3.7) 

JQ. 

Since thf‘ operator L' is with respect to = {x,y) and the above integral is w.r.to 77 , 
(3.71 can be written <is 

L'Poi:^ - C, A) + / L’Poiz - 77 , 77 , X)h{r), X)dr,n = 0 (3.8) 

d d ' 

Now consider L'Pq = [A — X^]Pq + a— + b— 4- c Pq 

where is the Dirac delta function and 


k{z,(A)= + + a(z-4.{.A) 


(3.9) 
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Now. (3.8) becomes. 

+ k{z.^. A) + + k{z,C, \)]h{r], X)d,,Q = 0 

->L (0 + k{z. A) + X)dr,Q + ^ k{z, e, A)]/i( 77, A)rf,0 = 0 

-dAO + k{z.^. X)-h{z,^, X) + J^kiz,^, X)]h(ri,^, X)dr,n = 0 


vSiiire f')A(0 = 0 for : 7^ we have 


h :,^,X) = k{z,^,X) + f k{z,r],X)hi7],^,X)dr,Q (3.10) 

JCl 

Ihe solution ot the integral equation (3.10) can be constructed by the successive ap- 
pro.Kimariun.s. Thus for /i(~,^,A) we have the series 

00 

h{z,^, A) = k{z, ■f, A) + k^{z, A) (3.11) 

i /=2 

where A:„{c,^,A) is the uth iteration of the kernel k{z,^,X). We show that the above 
.sorie.s is cunvcrgeiU w.r.to ^ in the region 0 for sufficiently large A. From equation(3.9), 


A-(c,TA) 


du ,du 
a — — f- b — — I- c 
ox ay 


1 /exp(— Ar)\ 

2^/2'kX \ ) 


where r -- |n - ^| = (x - + (y — ^ 2 Y- Then, we have 



cr + - ^i)a + (y - 6)^’] “ A[(x - ,fi)a + (y - ^ 2 )^] 


k { z ,^, X ) = 
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1 f exp(— Ar)^ 

27 s V ^ J 

”" (y ~ ^ 2 )^] - ^Vr[{x - ^ 1)0 + {y- 6)6] 


Using the continuity properties of a, 6 and c, we can find constants C[ and C 2 such 
that 


+ (2^) [(l^ - 6)0 + { y - 6)6) 
|\/f[(:c - 6)a + (y - 6)^] 


< C*; 

< C' 


From the above inequalities, it follows that 




1 exp(— Ar) 
\/A 5"2 


(C[ + AC^) 


1 exp(-fr) 
7A 7*2 


C{ exp 



+ AC 2 exp 


Fix Ao sufficiently large. Then for all A > Aq, we obtain 


exp 

|A:76A)| < C" 




or, 

exp (—^\z — 6) 

IM^,6A)| < C' -A {3.12) 

where C is a constant independent of A. (It can be noticed that for any positive r, 
[Aexp(— Ar)] — >■ 0 as A — >• 00 .) 


Now to obtain a bound for the 2nd iteration of the kernel k{z, 6 A), 
|A:2(z,6A)| = / k{z,T],X)kiT],^,X)dr,Q. 

J 0 



< 


(from (3.12)) 




drjQ 


< c exp(— (i.-ei)j^ 


< C" 


/2 


exp (-^\z-^\ 



exp (-rk-^l) 


n k-r7Pl77-eP 


dr,fi 


(3.13) 


To find a bound for the right hand side integral, we transform to polar co-ordinates, with 
origin at z and take one axis along the vector z — And letiz — if| = r-, \z — r]\ = ri 
and I?? — ^! = 7"2- 



Then, it follows that 


in |z - 77|2|77 - 


djjO. 


< 



exp (-Ari) 


sin 9 dO 
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= r dr, T' M 

Jo Jo — 2rri cos 9 


1 /‘^TT 

< -I exp (-Art) f 


sin 9 


1 + ^2 _ 2i cos 9 


d9 dt 


2 r°° 1 

= - exp(-Ari)— In 
r Jo ^ ^2t 


1 , 1 + + 2i 

l + t^-2i 


dt 


1 r°° 


r°° 1 1 4- f 

/ exp (-Art) - In dt 

Jo ^ ' t l-t 


< 


Ar2 


or. 


/p 




n ” A|z-eP 


(3.14) 


Using (3.14) in (3.13) we get 


C'2\ exp(^|z-{|) 



IM-.f,A)|<, |z-fP 

Now, we show by induction that for n > 2, 

|A„(z,<,A)|< ('yi p-§^exp(-^|z-{|^ 


(3.15) 


(3.16) 


• £ 

Assume that (3.16) is valid for n = m — 1 for some m, and we show that it is true or 


n = m. 
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\km{z,^, A)| 


/ km-i{z,rj,X)k{r},^,X)dr,Q. 
Ja 


■" k ix j ^ \z-r,P\v-e 



m — 2 


c 


’/2 


/ 

Jn 


ex p (xl^ - ^l) exp {- y\v - ^l) 

|z - 77|2|77 - .f|2 


drjQ. 


'C^Y V/2 f exp (^|z -r]\ + \\t] - ^l) exp (^|?7 - ^l) 

Ay Jn \z — vl'^lv — ^\‘^ 




< 



m~~2 


C'^exp(^|2-f| 


exp(^h-fl) 

/n \z — r)\‘^\ri — ^ 


From (3.14), we know that 

r !^pl±!ixz!)<i n < ^ 

J„ \z - n^v - (? " A|z-?p 

This leads to 

/^/2\’^-l fii2 f —X \ 
|A:m(2,e,A)| < \^—j jjr^exp |2:-,elj 


Hence the inequality (3.16). This inequality is valid for all A > Ao, where Aq is suffi- 
ciently large. Frona this the convergence of the (3.11) can be easily obtained. 
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3.3 Construction of the Fundamental Solution for L" 


In this section we construct the fundamental solution for the operator 


L" = 


. .odu , , 9 ,9 

A + aX^— + hX^— + cX^ 
ox ay 


or, L" = 


A + oA=|^+4A=|^i + {c+l)A= 

OX Oy 


X^ 


The procedure is same as in the previous section except that the kernel involved in 
the Fredholm integral equation is different in this case. The fundamental solution of 
the operator [A — A^] is given by (3.5). 

Following the lines of the previous section, we shall seek the fundamental solution 
in the form of (3.6). Then, the corresponding Fredholm integral equation is (3.10) but 
the kernel k{z, rj, A), in this case, is 

^o(^-e,e,A) 

The solution for h{z, f. A) in the integral equation (3.10) is same as (3.11) but the kernel 
and the corresponding iterated kernels will be given by (3.17). 

We show that the series (3.11) corresponding to the kernel (3.15) is convergent w.r.to 
^ in the region f2 for sufficiently large A. From (3.7) and (3.15) we have 


k{z,^,X) 


aA2|^^+6A"|^ + (c+l)A2 

OX ay 


k{z,^,X) = 


oA"|^ + i.A*|^ + (c+l)A2 
OX oy 


1 fexp{—Xr) 

2^[ , 
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where r — [z — — \J{x — (^i)^ + {y — ^ 2 )^- Then, we have 




1 /exp(— Ar)'\ 

27r\/A V / 

(c + l)A^r — ^A + — ^ {x — ^i)aA^ — ^A + — ^ {y — ^ 2 )bX^ 


_ 1 /exp(— Ar)A 

27r\/A V / 

|(c -r - 6)a -{y- 6)^]| - {( 2 ^ - 6)a + iv - 6)&} X^^/'l 

Using the continuit}’ properties of a, b and c, we can find constants C" and C 2 such 
that 


(c + l)ry‘+f ^ 


V2v/?. 


[(x - ^i)a + iy- ^ 2 )b] 


< C'l 


|\/^[(^ - 6)a + {y- 6)&]| < C'i 

Substituting these constants in the expression for we get 




r2 


1 exp 


(-td 


^A 


Cf A^ exp ( “ 2 ^ j ^^^2 exp 


A 


As in the previous section, for Aq sufficiently large, we get 


|A:(z,e,A)|<C" 


„ exp(-|| 2 r-.f|) 
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where C" is a constant independent of A. Further, in the similar way as in section-3.2, 
we can show that 

From this the convergence of the (3.9) can be easily obtained. 

3.4 Approximation of the Fundamental Solution 

The fundamental solution is given by the equation (3.6) where ^(77, (f, A) is represented 
by the convergent series (3.10). From (3.16) this series can be written by 

h{z, t A) = k(z, f, A) + x; k,(z,^. A) + O (a-') (3.17) 


for some suitable 1. Hence for large A, this can be written as 

h(z, f , A) ~ k[z, i,X) + 'h f. A) (3-18) 

j /=2 


The number of terms in this series could depends on A. For very large A, sometimes 
even k{z,^,X) will be a good approximation. To evaluate the right hand side integral 
in (3.6) we use a slightly modified version of the dual reciprocity method given in 
Partridge[54] as follows. The approximation for h{r],^,X) is proposed.by, t ' 

N+L 


(3.19) 
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where 4>3 are a set of initially unknown coefficients and the are approximation func- 
tions. N is the number of boundary nodes and L is the number of interior nodes in the 
problem domain. The above approximation is such that h is exact at these {N L) 
nodes. The expansion may then be considered as valid over the whole problem domain. 
Now, we choose a series of functions Uj such that 

[A - A2]% = /, (3.20) 


Now consider 

I Poiz-r],T],X)h{r},^,X)d^Q. = [ Pq(z-t],t],X) 

Ja Jci 


N+L 

J=1 


dr;^ 




(3.21) 


Since Pq is the fundamental solution of the operator [A - A^] we can use the boundary 
element formulation given in the section 2.4 of chapter 2. From (2.14), (3.21) can be 
written as 



Po h dO. 


N+L 


J=l 


^ r 

UijPYl / {dPQldn)u^dr 

k=i 
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Using the notations as given in section 2.4 in chapter 2, we get 

N-hL r N N 

HikUkj — GikQkj 
j=i L *=1 jfc=i 

N+L 

j=l 

This involves discretization of the boundary only. Internal nodes may be defined in the 
number and at the locations desired by the user; this is generally done at points where 
it is desirable to know the solutions. 

In order to find the vector 0, consider (3.19). By taking (N + L) different points, a set 
of equations like (3.19) is obtained. This may be expressed in matrix- vector form 

h = F(f) (3.22) 

where each column of F consists of a vector fj containing the values of function fj at 
{N + L) collocation points. As h is a known function, from (3.22) we have 

0 = F-'^h (3.23) 


Thus 4> will be a known vector. 
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3.5 Formulation of the method 

In this section, we show the formulation of the boundary element method method for 
the boundary value problem (3.1)-(3.3). The formulation for problem 2 will be similar. 
We integrate equation (3.1) with a weighting function w, 

(Lu)wdD. = J fwdO, (3.24) 



f {Au)wdQ, + [ a^wdCl + f b^wdCl 4- [ cuwdQ + f {—X^)uu:dQ. = [ fwdVl 
In Jn ox Jn ay Jn Jn Jn 


Using Green’s identities, 

dUfn 


[ (Aw)udQ + [ w^dT - [ u^dF - [ au^du - [ 

Jn Jr on Jr On Jn ox Jn 

dw 
dy 


da 


uwdFt 


+ f auwdy — [ hu^^dFl — f ^uwdQ — f 
Jr Jn ay Jn ay Jr 


n dx 

buwdx + / cuwdFl 
Jn 


X^uwdQ = j fwdQ 


L 


. dw .dw (da db \ 

Aw - a- b- \ ~ + — + c] w - X^w 

dx dy \dx dy J 

dw 


udCl + [ 

Jr dn 


[ u-:^dr + [ auwdy — [ buwdx — f fwdQ = 0 
Jr dn Jr Jr Jn 


[ (L*w)udQ- f ^udV+ [ awudy- [ bwudx- [ fwdQ = - / w^dF (3.25) 
Jn Jr dn Jr Jr Jn Jr dn 

where L* is the adjoint operator of L. 


i.e. L* 


A ^ 

^ ^dx ^dy 


'da db \ , 2 ' 
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Choosing w as the fundamental solution of the operator L* with respect to the point 
‘i’, where ‘z’ is any internal point of Q, where we want to find the solution. 

Then (3.25) becomes 

— u* — f ^^udT + [ awudy — f bwudx — f fwdQ, = — f w^^^dT (3.26) 

0Th JV JY J Cl JY OTt 

The expression (3.26) is valid for any internal point ‘i’ of Q, and it can be seen that the 
same expression is valid for a boundary point also. 

Now assume that is a polygonal domain and T = U"=i Tj- Further, u and (du/dn) 
are assumed to be constant on each Tj. Then for a boundary point we have 


-u* — [ ^^udr + f awudy — [ bwudx — f fwdQ = — f 

Jr on Jr Jr Jo Jr 


dUfji 


, w^dT 

r on 


-u‘ 


- [ ^^udT + [ awudy — f bwudx - f fwdQ = — / w^^ 

J[j^r, dn J[j^r, JU,rj Jo J[j^rj dn 


■dr 


-u 


y[ ^u^dT + Vf awUjdy-T^f bwujdx - f fwdQ = [ wq.dT 
yJrj on ■' jjr j 
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u 


^ [Jr, “ 4 + X X 


it;c?r 


u 


+ ^ 4- 5t — ^ Gy-gy 

J j 


(3.27) 


where = 


Gij 


X,I^‘*^+X,“”‘'^“X 


X, 


turfr 


bwdx 

k. Bi = f fwdQ. 
Jn 


Set ify = .ffy for z ^ 


1 + for i = j 


Then (3.27) gives 


n n 

yX HijUj -\- Bi = ^2 GijQj - (3.28) 

j j 

We know only some of the u^’s & g^’s. Keeping the known quantities on one side and 
the unknowns on the other we can write (3.28) in the matrix vector form 


Ax =■ g (3.29) 

Solving this system of linear equations we will know all the Uj's and g^-’s. Then again 
from (3.26) we have 

^ HijUj - Bi + J2 GijQj 

j 3 
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Thus we can find the solution at any internal point of Q, 


3.6 Numerical Examples 


Example 1 


. du , , 

Au + x— + y- X^u 

ox oy 


f , 0 < x,y < 1 & l<ICA<oo 


u = exp(— Arr) + exp(-Ay) on the boundary. 


Here / = -A[a:exp(— Aa:) + yexp(— Ay)]. This problem has the obvious exact solution. 
To evaluate Bi, we discretize the unit square into 16 square elements and then integ- 
rate / against the fundamental solution over each element using Gaussian quadrature 
formula and then take the sum of all these terms. One can avoid domain integration 
here by directly using the dual reciprocity method [4]. To evaluate the matrix elements 
Hij &: Gij we have used 4-point Gaussian quadrature formula. Gauss Elimination 
method (with partial pivoting) is used to solve the resulting linear system. 


To find out the solutions at the internal points, we define a grid of points {xk, yi) where 
Xk = k ! K V, k 1 ( 1 )^ yi = l/K + l;l = 1{1)K and JT > A so that we consider 
the solution at the points of the boundary layer region also. We call the approximate 
solution at the point {xk,yi) as Uki and define err{k,l) = — Uki\- Then the max- 

imum of this error function Max [err{k, Z)] which is tabulated for various A’s and 

k, I 

different number of boundary elements in table 3.1. As we use the asymptotic form for 
the fundamental solution we see that as A increases, the error decreases. And that we 
could achieve in just 4 elements! Also we note that the solution itself converges to zero 
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as A oo which coincides with the solution of the reduced problem. 


Table 3.1: Sup norm of the error function in Q of Probleml. 




No. of Boundary Elements 


A 

4 

8 

16 

32 

2 

1.952418 

1.203472 

.791541 

.483652 

10 

1.046732 

0.842572 

.197314 

.786261x10'^ 

20 

0.964321 

0.277501 

.960891x10-^ 

.342605x10-^ 

100 

.181534x10-^ 

.180921x10-^ 

.180701x10-^ 

.179072x10“^ 

200 

.921938x10-3 

.907138x10-3 

.894832x10-3 

.855028x10-3 

1000 

.628398x10-2° 

.593045x10-2° 

.563486x10-2° 

.552394x10-2° 

2000 

.732947x10-3° 

.709230x10-3° 

.690274x10-3° 

.690186x10-3° 


Example 2 




du du 
dx dy 




A 

dx 

dy 


“(0,y) = -M cos(7//2) 
u{x, 0) = —A, 


= /, 0<x,y<7r & l<cA<oo, 

— w(7r,y) = 0, 

and u{x, tt) = 0. 


Here / is given by 


/ = A^ exp(-A 7 r) - [M exp(-Mx) cos (y/2) + (1/2) exp (-Mx) sin (y/2) - A exp (-Ay)] ’ 
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where M = + (1/4). Then the exact solution is 

u{x,y) = exp (^-xsJX^ cos{y/2) + exp(-Ay) - exp(-A 7 r) 

Here the grid is taken as 2;fc = -Kk/K + l]k = l{l)K and yi = ttI/K+V, I = 1{\)K and 
K > X. The sup norms of the error function for various A’s are tabulated in table 3.2. 

Table 3.2: Sup norm of the error function in Q, of Problem2 


A 


No. of Boundary Elements 


4 

8 

16 

32 

2 

1.464029 

.983948 

.664352 

.329874 

10 

.963045 

.529384 

.279384 

.076261 

20 

.746531x10-^ 

.719387x10-^ 

.68390x10-^ 

.654158x10-^ 

100 

.453092x10-3 

.418373x10-3 

.384721x10-3 

.345127x10-3 

200 

.246721x10-3 

.228921x10-3 

.226018x10-3 

.222973x10-3 

1000 

.843278x10-2^ 

.837654x10-27 

.834396x10-27 

.831860x10-27 

2000 

.462931x10-^3 

.462745x10-^3 

.426691x10-^3 

.426672 xlO-'^® 


Example 3 


du du 

Au + Tr h — h U 

OX oy 


f , 0 < x,y < 1 & 0<eCl 


u — exp{—x/s) + exp{—y/s) on the boundary. 


Here f is given so as the problem has obvious exact solution. The numerical results are 
shown in table 3.3. 
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Table 3.3: Sup norm of the error function in Q of Problem 3 




No. of Boundary Elements 



4 

8 

16 

32 

10-1 

2.239343 

1.694729 

1.256905 

1.113904 

10-2 

.2418854x10-1 

.2418621x10-1 

.2418573x10-1 

.2417994x10-1 

10-3 

.4171974x10-2° 

.4171974x10-2° 

.4171974x10-2° 

.4171974x10-2° 

o 

1 

.311994x10-2°° 

.311994x10-2°° 

.311994x10-2°° 

.311994x10-2°° 


3.7 Conclusion 

In this chapter we have presented a boundary element method for general linear sin- 
gular perturbation problems of elliptic type. This we have done by constructing the 
fundamental solution corresponding to the given problem and have given an algorithm 
using the dual reciprocity method to evaluate the approximate fundamental solution 
through only boundary integrals. The numerical results show that the method has been 
successfully applied to various linear singular perturbation problems and the results are 


very encouraging. 



Chapter 4 

Iterative method for non-linear 
Singular Perturbation Problems 

4 . 1 Intro duction 

In this chapter, an iterative scheme is proposed for solving non- linear singular perturb- 
ation problem by linearizing the given problem. The linearized problem will then be 
solved by the approach discussed in the third chapter. The convergence of the iteration 
has also been established. The problem that we consider here will be of the form 

Lu = A.U = X^f{x,y,u) in 1 <C A < oo (4.1) 

u — u in r (= dQ) (4.2) 

We assume the conditions for the existence and uniqueness of the solution. That is, 
f(x, y, u) is a smooth function and there exists a constant /i > 0 such that fu{x, y, u) > 
/i-[53, 5, 33] 

4.2 Iterative method 

The quasilinearization can be done for the given equation (4.1) and the linearized prob- 
lems can then be solved using the approach discussed in the earlier chapters. But the 
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fundamental solution of the operator due to the quasilinearization happens to be de- 
pendent on the solution Um-i at the iteration, where Um-\ is the solution at the 
(m — 1)^'^ iteration. That is, at each iteration we have a different fundamental solu- 
tion. If we can have a common fundamental solution for all the iterations then that 
will save lot of computations. This motivates us to look for a scheme which leads us to 
have a fundamental solution independent of Um-i at the iteration. We propose the 
following linearization, and in the we prove that this iteration is convergent for large A. 

/, ju,,, = -Au,,! — X~Um ~ A [f 1 y ■) Um—i) n 7 ji_i/u( 2 ;, y, 'nn^— i)] in Cl (4.3) 

Uto = w on r (4-4) 

where wc' uo(^) as the initial guess solution and using this we iteratively solve for 

uiiz). Wi{z) etc. In all the iterations, the left hand side operator is [A - A^], 

and \\v know that the fundamental solution of the operator [A — A^] is 

p _ ^ exp(— A|z - ^1) 

° ^ 2v^ [ y'Alz -Tl J ^ 

lAsing this fundamental solution in the boundary element method, as discussed in 
chapter 3. we can solve (4.3)-(4.4). The formulation for this goes as follows. 

/ AumPo dCl- [ X'^UmPo dCl = f >?gm-m dCl (4.6) 

where y = [/(x,y,tx) - u/„<x,y,u)] and grn-i = [fi.^,y^Um-i)-Um-ifu{x,y,Um-i)]- 

We know from Green’s identity, 

j^£,u,„Poda = I^APaU„dn + ^ ^ 
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Using this in (4.6), we get 



.\s /)) is th(' funciamental solution of the operator [A — A^], we have [A — A^JPq = —^z, 
where A', is the Dirac delta function. Then, (4.8) becomes, 



The e<]uati()n (4.9) is valid for any internal point z of 0 and it can be seen that the 
same tfxpn'ssion is also valid for any boundary point. Hence for each boundary node 
4', we have 


/ dr = < + / dr+ [ x^gm-iPo dn ( 4 . 10 ) 

Jr dn Jr on Jo. 

We use till' dual reciprocity method[54] in order to transform the domain integral into 
buundary inlt'gral in the above equation. For this we first propose an approximation 
for as follows. 


fc., c. y (4.11) 

j=i 

where p^’s are a set of initially unknown coefficients and the fj are approximation 
functions. N is the number of boundary nodes and L is the number of interior nodes. 
The af)ove approximation is such that ^m-i is exact at these {N + L) nodes. The 
expansion may then be considered as valid over the whole problem domain. Now, we 

choose a series of functions % such that 
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[A - a2]z2, = /, 


(4.12) 


Substituting (4.11) in (4.10), we get 


Pc 


dlL-ni 




r dn - ' h dn Jn 


N+L 

E 

j=i 


dO, 




‘n+l 

i=i 


dQ 


r BP ' r 

+ I dr +11 ^ -^'^o[A - A2]u, dQ (4.13) 


\\> can again use the Green's identity and the property of the fundamental solution Pq, 
to transform above domain integrals into boundary integral. Therefore, after discretiz- 
ation (4.13) l)('Comes 



4" y~) Hij{Um) j 

3 


+ E 


Uij + 




3 


N+L , 

+ E 4’"“’ 


^ . iV 

Uij + ^ HikUkj — Y1 
k=l k=l 


N+L 


u: 




3=1 


(4.14) 


where 



dPo 


dn 


dr, 



L 


PodT, Qj 


duj 

HI' 




N 


N 




k-l 


k=l 
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(<hn)j = {du,Jdn)^, and for i^j 

= 1 + Hij for i = j 

Tht’ii (4.14) Ix'coines 

N+L 

^ ] j=l 

Nmt' tli.'it d.j is a known quantity. We shall explain, in a little while, how we find 
Sint'i' we ktiow on F. the right hand side of (4.15) is completely known. Then writing 
this tHjuntiou tor all the boundary nodes, we get a system in matrix-vector form. 

GQm = 6m (4.16) 

Hijlviiig this systtun, we can find all the (gm)j’s. We have seen that equation (4.10) is 
valid fur Ijoth int('nuil and boundary node. Hence equations (4.13) and (4.14) are also 
valid for any internal node ‘F. Writing the equation (4.14) for an internal node ‘r, 

- E 4 "''' \.«^i - 0 - 17 ) 

The right haiul side of (4.17) is completely known. Thus we can find the solution at 
any inttunal ntxift. 

la find th(! vector consider (4.11). Writing this in all the iV + L nodes, we get 

a .system in mat rix-v(>ctor form 




(4.18) 
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where each column of F consists of a vector fj containing the values of function fj at 
{N + L) collocation points. As Qm-i is a known function, from (3.22) we have 


^ (4.19) 

Thus can be found. This procedure involves discretization of the boundary 

only. The number and the locations of the internal nodes may be defined as desired; 
this is generally done at the points where it is desirable to know the solutions. 


4.3 Convergence 

To see the convergence of the iterative method, we consider the equation (4.9). That is, 
for z G we have 

Umiz) = - ^ X^9m-iPo dQ + Po^ dT - ^Um dV (4.20) 
For exact solution u{z) we have from equation (4.1), 


Au — \^u = A^(/ — u) 


And integrating the above equation against the fundamental solution Pq, 


JjAu - X\]dQ = X\f - u)Podn 
From (4.6), (4.7), (4.8) and (4.9) we get 

u(.) = - V(/ - u)Po dn + dr - <jr 


Subtracting the exact solution u from Um, we get 

dum du 
dn dn 


Um - U 


A^Fobm-i - (/ - w)] dQ + J^Po 


dT 

(4.21) 



55 


Since on F, u =z Um, the third integral on the right hand side of (4.14) goes to zero. 
This leads (4.14) to 


Um-u = 



if -u)]dn + J^P 


dUm 

dn 


du 

dn 


dT 


Um-U = 


OiJm-l 


U 


m— 1 


ifu)r 


f + u)]dn + j^P 


dUm 

dn 


du 

dn 


dr 



f '^m—lifu')m—l d” '^fu 


ufu + ■w)] dQ + J^P 


dujji 

dn 


du 

dn 


dV 


The Mean Value Theorem gives 


fm—l f — /u(fi) i'^m—1 '^) S'Ud [nm—lfui^m—l) ^/u(^)] — ['^fuui^')~^ fu]i'^}i'^m—l ”^) 

where u = tUm-i + (1 — t)'^, 0 <t <l 


Choose Mm such that 


|/u({t)| < Mm 


and 

lfi/uu(fi) d“ yu(fi)| — Mm 


Then we get, 


l/m-l - / - ^im-l/«(^^m-l) + ufu{u)\ = ^^m \Um-l - u\ 
Also, let Lm be a constant such that 

\ufuiu) -u\<Lm 


So, (4.21) becomes, 

< J^X'^PQiz-^,f,X)[2Mm\Um-l-u\+Lm] dn + J^Pq 


\Um - u 


dUm du 


dn dn 


dT (4.22) 
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Consider the second integral on the right hand side of (4.22). Since it is over the 
boundary T, we have | 2 : - ,f| > 0. Hence for very large A, Po(^ - A) will be very 
small. That implies that the integral value becomes very very small, i.e., 

I 


f A) 

Jv 


dn 


du 

dn 


dT < s', for some small s' 


(4.23) 


Now let us consider the first integral term. 


/ A^[2A/m|Um-l — u| + Lm]Po dQ^ 

•J iQ 


+ 


[ X‘^[2M^\Urn 

jQ\Cp 

[ X^[2Mrr,\Um-l 


“1 4“ Ijxx}\P(i 

lil + TjtiJTo dfl^4.24) 


where Cp is a circle of radius p with centre z for some small p > 0. Again, in the first 
integral on the right hand side of (4.24), |z — ^| > 0 which means that value of this 
integral is very small for A very large. We have therefore. 


/ )?PQ\2Mm |unj-i — u| + Lr„\ dQf < s" for some small s" 
Jn\Cp 

For integral over the circle Cp, 


(4.25) 


A' 


■Po[2Mm\Um-l-u\+Lm] dO.^ = X‘^[2Mm \Um-l{z) - u{z)\+Lm] / Pq{z-^,C, X)d^ 


(4.26) 


I = i^. 


2V^ 


exp(-A|z - ^\) 


1 c- 

/%r Jo Jo 


2v^ 


r exp(- 


■Xt 


L \/^ 


d^ 

t dt 


Furthermore, 
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where we have taken the co-ordinate transformation (<fi -x) = t cos 9, (^2 -y) =t sin 6. 


fpo(z-u.m = ^r\?^]tdt 

Op 2\/ 27r Jo yXt 


2\/X Jo 


cp(—Xt)Vi dt 


By mean value theorem for integrals, 


P(){z — X'd^ = ^^'^pexp{—Xp')yf^ where 0 < p' < p. 


Substituting this value in (4.24), we get 


A^Po [2Mm\um-iu\ -r Lm] dQ^ = p exp{Xp')^ \um-i{z) - u\+X^/^ p exp{Xp') Lr. 


For given A, p can be chosen such that. 


< r < 1 & \Lr, 


for some s'" very small. Then (4.26) becomes. 


0^/n' < 5 '" 


(4.27) 


(4.28) 


A'Po [ 2 Mm 


Um-l - u| -f- Lm] dQ^ < r \Um-liz) - u| + s'" (4.29) 


This leads (4.24) to 


f Mm>^^\Um-l - u\Po dO,^ < s" -f s'" + T \Um-l ” w| 
Jo 


(4.30) 


Now, from (4.22), (4.23) and (4.30) we have 


Um-nl < s' -1- s" -f s'" ■+- r \Um-l - U 


u| < s 4- r \um-i - u\ where s = s' + s" 4- s'" (4-31) 
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Then by induction on m, we have 

\um - u| < r”"|uo - u| + (1 + r + + . . . + (4.32) 

which ensures the convergence of {^^( 2 ^)}- 

4.4 Test Problem 

As an example we solve the following problem [6]. 

Aw = \^[u — g^u^] in Q. : 0 < a:, j/ < tt; 1 ^ A < 00 . 

14 = 0 in 50 

The function g{x, y) is determined so as the exact solution is 

u(x, y) = a;[exp(— Aa;) — exp(— A7r]sm(y) 

Here the reduced problem, that is the problem obtained by setting 1/A = 0, has the 
solution u = 1/g. The linearized problems are 

Aujn - = A^ in 0; 1 < A < 00 

Urn = 0 in 50 


We can solve the above problem and find the solution at any internal point. In this 
example we solve the problem by using 4, 8, 16 boundary elements. In the domain O to 
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Table 4.1; Sup norm of the error function in Q 


No. of Boundary Elements 
4 8 16 


10 

.965 

(21) 

.839 

(19) 

.673 

(18) 

20 

.450 

(19) 

.293 

(18) 

.209 

(18) 

100 

.964x10“^ 

(15) 

.737x10-1 

(14) 

.605x10-1 

(14) 

200 

.405x10-1 

(15) 

.318x10-1 

(15) 

.247x10-1 

(14) 

1000 

.647x10-3 

(9) 

.647x10-3 

(9) 

.647x10-3 

(19) 

2000 

.483x10-^ 

(6) 

.483x10-'^ 

(6) 

.483x10-^ 

(6) 
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the fact that we have actually constructed the fundamental solution for general linear 
olli ptic problems in chapter 3. 



Chapter 5 

Application of the method in 
Semiconductor Physics 

5 . 1 Introduction 

In this chapter, we propose to apply our technique to a problem in semi-conductor 
physics. In earlier chapters, we have seen how our method can be applied to both 
linear and non-linear problems. Here, we solve a non-linear system of partial differential 
equations using the iterative technique proposed in chapter 4. 

The fundamental semiconductor device equations form a system of three second 
order elliptic differential equations subject to mixed Neumann-Dirichlet boundary con- 
ditions. The system consists of Poisson’s equation and the continuity equations and 
describes potential and carrier distribution in an arbitrary semiconductor device. 

Let i) denotes the electrostatic potential, n denotes the electron density and p de- 
notes the hole density. Then the basic semiconductor device equations are[48] 


SsAi; = q(n-p- C{x,y)) 

(5.1) 

div(D„Vn - /.i„nVV>) = 0 

(5.2) 

div(DpVp -f pppVtl}) = 0 

(5.3) 
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{x, y) € il, where ia a bounded domain in representing the device geometry. The 
equation (5.1) is the Poisson’s equation and equations (5.2) and (5.3) are the electron 
continuity equation and hole continuity equation respectively. The other variables in 
the problem are 


£5 

; semiconductor permittivity 

Q 

elementary charge 

Dn 

electron diffusion coefficient 

Dp 

hole diffusion coefficient 

/in 

electron mobility 

flp 

: hole mobility 

C{x, y) : 

doping profile (i.e., the difference of the electrically active concentration 

of donors and the electrically active concentration of acceptors.) 


The boundary dQ splits up into three adjacent parts, namely and dQos 

such that = U3t=s3 Ck where Ck are connected arcs and CiflCj = { } for i ^ j. 
Dirichlet boundarv' conditions for n, p are given on dClc as follows. 

{n-p- C{x, y)) lane = 0 and npjane = nf 

where is the intrinsic number of the semiconductor. For the potential rp, we have 

7P\c,=UTln- + Uk 

Til 

where Uk represents the potential applied to the ohmic contacts Ck, and Ut is the 
thermal voltage. 
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5.2 Scaling of the Problem 

By appropriate scaling, the above problem can be viewed as a singularly perturbed one 
in the following way. Assume the validity of Einstein’s relations 

■^n jj 

l-^n t^p 

And take the transformation called Boltzmann statistics, 

n = n,- exp{'ip/UT)u, p = nj exp(— 

where u = exp(— (p„/17r), v = exg{4)p/U'r) and (j)p are the electron and hole 
quasifermilevels respectively and also u > 0, u > 0. Then, (5.1) becomes, 

= g(n, exp(^/£/T)w - Ui exp(-^/17r)u - C{x, y)) (5.4) 


Again, from Boltzmann statistics, we have, 


DnVn - yn-nVip = Uryn^irii exp(V’/C/r)u) - y-nrii exp{Tp/UT)uVip 


= Ury-n'i^i^ /U t)u) - ij.nniexp{ip/UT)uV‘ip 


Vib 

u ex'p{^ IU t)— — 1- exp(^/C/j')Vu 
Ut 


= Ury-nrii 
= C/TPnniexp(^/C/T)Vu 


— pn, exp{'ip/UT)uV'ip 


Therefore, from (5.2) we have, 


div(Z)„Vn — fXnTiV'ip) = 0 
div(t/Vp„niexp('i/)/i7T)Vu) = 0 
div(pnexp('0/17r)Vu) = 0 


(5.5) 



Following the similar lines, we have 


DpVp + PppV-ip = 


= UTPpV{ni exp(-^//7r)n) + ppm exp(-V>/C/r)t;V^ 
= UtPpUiV ( exp{—il;/UT)v) + Ppn^ exp(-'ip/UT)v'V'ip 

U'J' PpTt^i 

+ ppniexp(—'ip/UT)vV'ip 
= UTMpniexp(-il^/UT)'Vv 


—V exp{-t]j /U t)^^ + exp{-'tp/UT)'Vv 

Ux 


And, hence from (5.3) we have, 

div(Dp Vp 4- Ppp'V'ip) = 0 
div(/7rPp7ii exp(— •0/177’) Vn) = 0 
div(/ipexp(— 0 / 1 / 7 ’) Vu) = 0 


Assume that C is bounded in Q and set 

C 

C = sup |C'(a:, y)|, D = and Z = diam(Sl). 

n C 

The dependent variables are scaled as follows. 

0 n p 

And the independent variables are scaled as 

X y r / \ A 

3^5 = y, for (x^.y^) 6 iZ 

Then, equations (5.3)-(5.5) are transformed to 

e^A0 = (i^[exp(0)'u — exp(— 0)u] — D 
div(exp(0)Vu) = 0 
div(exp(-0)Vt;) = 0 
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or, 


e^A'0 = 5^[exp(^)w 
d‘6 du dip du 
dxdx^ dy dy 


Av 


dip dv dip dv 
dx dx dy dy 


exp(~ip)v] — D 
0 

0 


( 5 . 7 ) 

( 5 . 8 ) 

( 5 . 9 ) 


where 


/ J - PqC' 


and — 


Til 

0 


Here is the minimal Debye length of the device. The scaled boundary conditions are 


dip 

du 

dv 

dn 

«n„ ” ^ 

an. ” 

dip 

du 

dv 

dn 




(5.10) 


( 5 . 11 ) 


hIc. = e.xp(-(7i/Ur), v|c. = exp (!7i/!7r). (5.12) 


V'Ict 


In 


252 


Ck 


Vt 


(5.13) 


For modern devices 0 > 10^^ cm^. For silicon and silicon oxide at room temperature 
T = 300i^, we have q = 10"^® As, e. = lO-^^^s/Vcm, Ut = 0.025F. And with real- 
istic value / = 5 X cm, we get s < 10"^ < 1. Therefore, the problem (5.7)-(5.13) 
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constitutes a singularly perturbed quasilinear elliptic system of differential equations. 


Normally, the parameter < 10 ^ <C 1, too. This however gets compensated by 
the Dirichlet boundary conditions (5.12) and (5.13) which imply that, on Ck, 


(5^(exp('^)n - exp(— '0)r) = < exp [ In 


exp ( — In 


D + VD^ + 454 


2(J2 




D + VW+W 

252 


El 

Ut^ 


V 



Uk 




Uk 


Uk 


exp I --y exp 


Uk 


d + vd^TW 


252 


d + ^/d^TW 




2^2 


D + \/02 + 454 


25^ 

2 


D + ^/D^■¥A5^_ 


= 0(1) as 5^ — )• 0 

Hence, we regard as a fixed parameter and solve the problem (5.7)-(5.13) keeping s 
as the perturbation parameter, i.e., letting e — > 0+. 


5.3 The Iterative Method 

In this section, we propose an iterative technique, an extension of the one we have used 
in chapter 4, to solve the singularly perturbed quasilinear elliptic system (5.7)-(5.13) 


Integrate the equation (5.7), with a weighting function wi, we get 
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f dQ= f [S'^exp{ip)u - 5^exp{-'ip)v - D]wi dO. 

*'0 «/ o 

Set A = 1/e, then we have, 

/ dQ, = f X^[S^ exp{'ip)u — (5^exp(— •0)t; — D]wi dO, 

Cl J Cl 

Subtracting / X'^'ipwi dQ on both sides, 

Jci 


(5.14) 


[ [Alp — X^-tpjwi d9.= f exp('0)u — exp(— — D — ip]wi dO. (5.15) 
Jn Jn 

From Green’s identity, we have 

J Alp Wi dfl = j Awi Ip dQ + J ^ ~ J where F = 5r2 (5.16) 

Using this in (5.15), we get 

[ [Arui — A^md v dCl + [ dT — [ tp-^ dT 

Jq Jr on Jr on 

= f X^[6~ exp{'tp)u — 5^ exp{—'tp)v — D — tp]wi d9t 
Jo. 


Let z = {x,y) and ^ = (^j'fa) be points on Q. Let us choose Wi as the fundamental 
solution of the operator [A — A^], then, 

= / A^[d^ exp(’0)u - 5^ exp(-'0)u — D — 'ipjwi dO. 

Jn 

(where J^(^) is the Dirac delta function) 
Then, by the property of Dirac delta function, 

9v ... _ /v-^dT 

Jr on 

+ f A^[(5^ exp('0)u — 6^ exp{—'ip)v — D — 'ip]wi dfi(5.17) 

Jn 
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This is the representation we get for V’(-z). For u{z), integrate equation (5.8) with a 
weighting function 102 , 


j (Au)w 2 dQ + J 


dip du 
n dx dx 


f dipdu 

wdQ+ —~w dfl 
Jo. dy dy 


From Green’s identity (5.16), 

du „ f dw 2 f dip dw 2 


[(Au„)uda+fw,pdr-fu^dr+[^u?p-dci - 

Ja Jr dn Jr dn Jn ox dx Jq ox^ 


+ f ^UW 2 dy- f 
Jr dx Jci 


dip dw2 


u- 


n dy dy 


dO. 


-L 


d'^ip 
h dy'^ 


UW 2 dO, 


-I 


n dx^ 
dip 


T dy 


UW 2 dx = 0 


f . dipdw2 dipdw2 . jr, , f JT 

/ AW2 - 5 ^ ^ AvW2 U dQ + / W2^ dr 

Jn dx dx dy dy Jr on 


dy dy 

f dw2 , f dip 
- / u-^ dr+ — uu ;2 dy 

Jr dn Jr dx 


-I 


dip 
r dy 


UW 2 dx = 0 


Choosing W 2 as the fundamental solution of the operator 


dip d 
dx dx 


dip d 
dy dy 


Alp 


we get, 

, , f dw2 ^ f dip j ^ f dip , fdu 

or, 


u{z) 




Similarly, we have for v{z), 
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where wz is the fundamental solution of the operator 


A + 


di}) d 
dx dx 


dtp d 
dy dy 


+ 


Let us assume, for the time being, that we know tp, (dtp/dn), u, (du/dn), v , , (dv/dn) 
on whole of F. And let tpo(z) be the initial approximation for tp(z). Then from (5.18) 
and (5.19) we have. 




and, 


^^1 


= - /r - L ^ Jr /r "“i 


Then from (5.17), 

tP,iz) = (uie^° -uxe-’^°)-D-V'oKdF! (5.22) 


We thus define the iterations. 


u. 


.w = -/. 




r dn 


u 


m— 1) 


/r <‘^+1 an 

(5.23) 


dr 


m~l) 


= - /r If ” Wr a„ 

(5.24) 


(ir 



71 


Wm{^) = / (^^me’*'”“' - Ume~^’"~')-i^-V’m-l]iyi 

(5.25) 

We know from cjhapter-4 that the sequence {^m(.z)}“=i converges to tp{z), the solution 
of (5.7). 

5.4 Boundary Element Discretization 

To find the solution ipiz), for .z € we first define N nodes on the boundary and 
discretize the boundary into N elements, i.e. T = U^i Ty. Then equations (5.23)- 
(5.25) becomes, 




r 8 u .2 

Ur..r, Sn 


T-p f (m— 1) 1 


dx+f 

".r, 3y Jul.r, 


(m-1) du 


M = -[ 




Uf=ir, dn 


-i 


f d'lprn—l (tti— 1) 7 

/■„ + 

'U",r, 9!/ 


'ur^r, 


^m{z) 



H - '!/),n_i]iyi dO 
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or, 


Um{z) 


^ r 


-7r, an a^r 

du 
dn 




u dy 


dV 


(5.26) 




(5.27) 




Call 



+ A^[(5^ (ume^'^~'- - Vme - D - tpm-llwi dCl 


(5.28) 


6 ^_l = \^[5^ - Vrae-'^-A -D- i>m~l]w, dQ 

which can be evaluated through boundary integrals using the dual reciprocity method 
as discussed in chapier-3. Now, to evaluate integrals over T^, we can approximate Fj 
to be linear, quadratic or cubic elements. For instance, if we take it to be quadratic, 
then on each Fj, we have 3 nodes, say, {xi,yi), {x 2 ,y 2 ) and (^ 3 ,^ 3 ). And the functions 
are approximated on Fj as follows. 


ip{s) 

dtpis) 

dn 

u{s) 

du{s) 

dn 



(5.29) 

(5.30) 

(5.31) 

(5.32) 
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where 


u(s) 

dv(s) 

dn 





C=1 



(5.33) 

(5.34) 


Niis) = -^(l-s), A^2(5) = (l+5)(l-s) and A/ 3 (s) = |(H-s); -1 < s < 1. 

And i^cj is the value of ^ at the node c of the element and (d'ip/dn)c,j is the value 
of (dTp/dn) at the node c of F^ and so on. The spatial co-ordinates are related by the 
transformation, 


^(■®) = Ncis)xc and y{s) - ^ Nc{s)yc 


C=l 


C=1 


In the boundary integrals, 




where 




Now, substituting equations (5.29)-(5.34) in equations (5.26)-(5.28), we get 
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Rearranging the summation terms, we get 


N 3 


Um{^) — f 

1 1 X. 




j=lc=l 
N 3 


r, dn 


^ Nc{s) dT 


+ 


EE«o, L ^wr-^^N.(s) dy 

j-lc-l 


j=l C=1 

N 3 


+ 


SS S./r 


wt~^'>Nc{s) dT 


Vmiz) 


N 3 f. 

EE-c, L 

j^l c=l -'ry 


dn 


Nc{s) dT 


( 5 . 35 ) 


( 5 . 36 ) 


( 5 . 37 ) 


( 5 . 38 ) 
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+ 


j=ic=i ^y 




^mi^) — 


* SSliJJ, 


£r(g) / ..«.} dr 

J=:lC=l \^'^y C,J 


N 3 




j=l C=1 


(5.39) 


(5.40) 


Set 




( 2 ) 


y>c 


/r, It + 1 ^»r-‘'iVe(.) * - 1 .. 




(2) 


ij,c 


= 


Tj (9n 

^ ?/;f-'’A;(s) dr 
<9u-, 




‘ZJ,C 


1 . (m— 1) AT / \ T . /" ^'^Tn—1 (m— 1) f \ j 

^ iVc(s) dar 


g: 


( 3 ) 

ij,c 


^(I) 

U>c 

r-(l) 

^d.c 


- / ^ iv,(5) dr - / f 

Jrj dn Jr, dx ^ ^ ^ dr 

^ ti;f-'d\;(5) dr 

^ u;iA^,( 5! dr 


r, 5y 


where the subscript ‘f denotes the internal node 2 : where we find the solution. Note that 

the right hand sides of the above expressions are all known, hence all the h\]\, 

* 

^i/,ci 3.re known entities. Now, the equations (5.38)-(5.40) become. 


u. 


= 


N 3 

EE*-. 

J = lC=:l 


«<l+EE 

j=l c=l 





( 5 . 41 ) 
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Vr 




N 3 

El 

7 = 1 c=:l 

N 3 


.W = EE»«^S + EEf^) Gg 

■ ■ j=lc=l c,j 


w = ee(I!) <>-ee^.xu»^-. 


j=ic=i \dn) 


C,J 


J=zl C =1 


(5.42) 


(5.43) 


On Fj, we have (Fig-5.1) 


r— { 

1 

II 

t-F 

U2J — Ujj 

^3,i 

— Ujj-i 


II 

^3,j 

= Uj+l 


V’2,i = Uj, 


= ^J+l J 


Similarly, 


'§3l\ 






\dnjj > \dnj^j \dnjj^i 


( f dv'\ ( d2i\ ( dri\ V 

\dnJi^ VSn/j-l ’ \ 9 njj’ {drij^j \dnjj^i 


fM) =fM) (M) —(§±] (§t) z=:(^) 

\3nyy_i’ \d 71 J 2 j \dn)j'' \dnJ3j \dn)j^i ^ 


(5.44) 


(5.45) 


Thus we can find Um{z),Vm{z), and i^miz) for any internal node z 6 provided 
we know ail the Uj,{du/dn)j, Vj,{dv/dn)j, ipj,{dijj/dn)j. But, note that we know 
{dipfdn), (du/dn) and (dv/dn) only on an.d we know u,v and only 

on UfcOfc. That is, we do not know {du/dn)j, {dv/dn)j, and {d'ijj/dn)j for some j’s 
and ipj, Uj and Vj for the remaining y’s. Now, to find these unknown quantities, we go 



Figure 5.1: The (boundary) element Fj 

back to the equations (5.23)-(5.25). We can see that these equations are valid for any 
boundary point also. Thus, equations (5.26)-(5.28) and equations (5.41)-(5.43) axe also 
valid for any boundary node i. From (5.44) and (5.45), equation (5.41)-(5.43) can be 
written as 


■Wo, [w(2) . 0-(2) , rr{2)] \ [^(2) , ^(2) p(2) 

2^ Uj + iaiy,3j + ^ j + ^ij,2 -b ^ij,z 


(5.46) 


i; V, [jjgl + + E ( . NS + Ggi + G^: 


( 5 . 47 ) 
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A = Efif) [GS+GaUG<;y-EV'4<i+jygUH™]+6„..(5.48) 

j=i j j=i 


H<;> = [Hgi + 4‘i + 4‘i] 
4 > = [ 41 + 41 + 41 ] 
4 > = [ 41 + 41 + 41 ] 


Gil* - [gSIi + gIIj + 4,3 

G144+41+GS] 

r^(3) _ rr^(3) , ^(3) , ^(3) ] 
^ij — I + ^ij,2 + 


Then writing the equations (5.46)-(5.48) for all the boundary nodes, 


E<4 = E 

j=i j=i 


/^(2) 

h\^L 

j=:i \ / m 


i-E<«. 

3 = 1 


j ry{3) _ 


^=1 / m 


i=i 


= Ef^VG'l* + i.-. 




= l-H 


1 - 


Hf = 1 + 4'’ 


for i j 




Hlf = 4 > 


for i 7 ^ j 


Then, we have 


E<«S’ 

j-1 


^ Gif 

£<[dnj^ 

^ = 1 \ / 771 


(5.49) 


i=i 


y f g^3) 


(5.50) 


E44 


E f^T Gll> 

■Adn)^ ^ 

j=l \ / m 


4' + i>m-l 


(5.51) 
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Keeping the known quantities on one side and the unknowns on the other, we get three 
systems of equations, which can be symbolically represented by 

1 



(5.52) 


= h(3) 

m m m J 

Solving these linear systems, we can find all the 
{dtpi^/dny^ for all j. 


5.5 Numerical Results 


We have solved a two dimensional diode and presented here the numerical results. The 
geometry of the device is 

Q = {(x,y) : 0 < x,y < 1}; flo = + = n\Oo 


C = 


10^^ cm ^ in fii 
-10^^ cm~^ in fio 


; = 5 X 10 ^ cm 


The boundary dQc is splitted as 

Cl = |(a:, 0):0<a:<i|; €2 = {( 2 :, 1) : 0 < a; < 1} 


Here = 52 _ j^q-t problem[23] as discussed in the section 5.3 and 

5.4. Fig-5.3 shows the potential in thermal equilibrium; i.e. Uq = Ui = OV. Fig-5.4 
shows the electron density n and Fig-5.5 shows the hole density p. These results have 
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Figure 5.2: Geometry of the 2-D diode 

been obtained in 20 iterations and using just 32 boundary elements. 

The entities which are defined in terms of 

the boundary integrals are evaluated using 4-point Gaussian quadrature formula and 
resultant matrix systems (5.52) have been solved using Gauss Elimination method (with 
partial pivoting). 


5.6 Conclusion 

As it stands, boundary element method has very many advantages over domain dis- 
cretization methods as it discretize only the boundary rather than domain. And, since 
we have proposed an iterative boundary element method which uses a fundamental 
solution in the asymptotic form, it efficiently solves the singularly perturbed problems. 
Further, the method is simple and straightforward and does not require the know- 
ledge of the behavior of the solution expected, hence is very much suitable for those 
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non-mathematicians who do not appreciate the use of mathematical tools to obtain the 
qualitative behavior of the solution system. 
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